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p ' Abstract 



A numerical analysis of mechanical frictionless contact between rough self- 
affine elastic manifolds was carried out. It is shown that the lower cutoff 
wavenumber in surface spectra is a key parameter controlling the represen- 
tativity of the numerical model. Using this notion we demonstrate that for 
representative surfaces the evolution of the real contact area with load is 
universal and independent of the Hurst roughness exponent. By introduc- 
ing a universal law containing three constants, we extend the study of this 
evolution beyond the limit of infinitesimal area fractions. 



1 Introduction 

Real surfaces are self-affine [H [2] and the surface heights are distributed nor- 
mally [5]. Due to this roughness the real contact area A is often only a 
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small fraction of the nominal or apparent contact area Aq. The real con- 
tact area fraction A/Aq and its evolution determines the contact resistivity, 
friction and the transfer of energy (heat and electric charge) through the con- 
tact interface. Thus understanding the contact-area evolution has profound 
implications in various fundamental (e.g., origin of friction) and engineer- 
ing studies (e.g., electro-mechanical contact, tire-road interaction). From 
early experiments [1] and analytical theories [31 E], it has been considered 
as an established knowledge that A evolves linearly with the normal load 
F at relatively small fractions of contact. It was confirmed by numerical 
simulations of normal frictionless elastic contact between rough surfaces [6], 
which made considerable progress since then. Nowadays, the most advanced 
computations exploit fine surface discretizations (2048 x 2048 [3 [HI 13 ITU] . 
4096 X 4096 [11]) to extract statistically meaningful results valid even for 
small fractions of contact. Many of the numerical investigations are based 
on synthesized surfaces preserving the aspect of self-affinity. In early studies 
[6], [12] the artificial rough surfaces obeyed fractality which scaled down to 
the domain discretization [13]. Later it was realized that the surface has to 
be smooth enough [TJ [T3] to represent correctly the mechanics of contact and 
to obtain a reliable estimation of the contact area growth comparable with 
analytical theories. This smoothness implies a sharp decay of the power spec- 
trum density for wavelengths smaller than short cutoff wavelength A^. On 
the other hand, to the best of our knowledge, there are no consistent study 
on the influence of long cutoff wavelength A; on the response of rough sur- 
faces. Despite remarks made in [15] on the importance of A; for mechanics of 
rough contact, in almost all investigations of the real-contact-area evolution 
[3 [To], [16] A; was limited by the size of the specimen A; = L. In this article 
we show that A; plays a crucial role in numerically precise theory of rough 
contact, because it controls the representativity of the simulated system. 

In the limit of inflnitesimal contact, analytical theories predict that the 
real contact area A evolves proportionally to the normal force F with coeffi- 
cient K normalized by the root mean squared slope of the surface \/J}Vh\^ 
and the effective Young's modulus E* [T7], which gives 

A = k/./{\vW)F/E*. (1) 

According to asperity-based [3] models (Bush, Gibson, Thomas (BGT) [5], 
Greenwood [18], Nayak-Thomas [19]) k = \/2it but Eq. ([1]) is valid only in 
an asymptotic limit for infinitesimal fractions of the contact area A/Aq — )■ 0. 
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FIG. 1: Representativity of a rough self-affine surface can be estimated accord- 
ing to the normahty of heights distribution. Here, the criterion of proximity is 
an error between the distribution of surface's heights and the associated normal 
distribution; (a) this error is strongly dependent on the cutoff wavenumbers kg 
and ki. Every point on the graph represents an error averaged over 30 statis- 
tically equivalent surfaces. Examples of rough surfaces and the corresponding 
heights distributions for different cutoffs ks = 40 and (b) ki = 1, (c) ki = 4, (d) 
ki = 16 are also depicted. 



3 



The convergence of k, with decreasing fraction A/Aq is slow and depends 
on the bandwidth parameter a introduced in [20] as a = mQm4^/m\, where 
rrii = J^' k^C{k)dk is the z-th moment of the power spectrum density C{k) 
of the surface. For instance, for a = 10, k is overestimated by about 12% 
for A/Aq = 10~^ and this error reduces only to 4% when A/Aq = 10~^° [18] . 
Consequently, using a complete numerical model it is not possible to demon- 
strate the linearity predicted by Eq. ([T]) nor to approach its asymptotic limit. 
It is also worth noting that in forementioned models the evolution of the real 
contact area is strictly nonlinear for realistic fractions of contact |2T1 122] . 
Moreover, this evolution does not depend on the Hurst roughness exponent 
H but only on a |2l] . A competing theory was proposed by Persson [231 121] , 
it also predicts a linear contact area evolution for small contact fractions, 
however, the obtained coefficient of proportionality k = ^/sJtt is signifi- 
cantly smaller than in the asperity-based models. As in the latter models, 
the evolution of the contact area does not depend on the Hurst exponent. 
Numerical studies [H [TJ [TOl [16] demonstrated an approximately linear evolu- 
tion of the contact area with load. But in contrast to the analytical models, k 
was shown to depend on the Hurst exponent and to be confined between the 
asymptotic limits of the BGT [5] and Persson [231 IM] theories. In this article, 
by introducing the notion of surface representativity, we obtain qualitatively 
new numerical results, which we believe correctly represent the mechanical 
response of realistic rough surfaces. Moreover, they appeal to a broad in- 
terpretation of the data not restricted to a single proportionality coefficient 
between contact area and load for infinitesimal contact fractions. 

2 Representativity of rough surfaces 

A meaningful numerical simulation of rough contact has to be carried out 
on a representative self-affine surface element (RSSE) (e.g. [25]) either gen- 
erated numerically or chosen from experimental measurements. The reason 
to use generated rough surfaces is that they may be obtained for any cutoff 
wavelengths and Hurst exponent H. On the one hand, the RSSE has to be 
large enough to obtain a similar response for different realizations of statisti- 
cally equivalent surfaces. On the other hand, it has to be as small as possible 
to retain a numerically solvable contact problem. So the RSSE is defined 
according to the permitted error between mechanical responses of surfaces 
of different sizes, e.g. L and 2L. This mechanical representativity can be 



4 



12 



10 



c 

o 



u 

O 

B 4 
c 
o 
o 



Cutoffs 
L/A; L/A, 
I- 16 32 
16 128 
1 32 
1 128 



Hurst exponent 
0.20 0.36 0.52 0.68 0.84 

i i i f i 

□ o S ^ 5 

S 5 § 5 5 ^ 
S i * f 5 ^1 



1 .i 









4lh¥ 



Asymptotics, BGT [5] 
Contact evolution law, Eq. (2) 
Greenwood [18] 
BGT [5, 21] 

Asymptotics, Persson [23] 



0.0 0.01 0.02 0.03 0.04 0.05 

Normalized external pressure, pg/£* 



FIG. 2: Evolution of the real contact area fraction A/Aq with normalized ex- 
ternal load po/E* are plotted for different cutoff wavelengths and Hurst expo- 
nents; the rms slope is kept constant for all surfaces (|V/iP) = 0.1. The re- 
sults for representative surfaces (A; = L/\Q) are in a good agreement with the 
asymptotic limit of the Bush, Gibson, Thomas (BGT) theory [5] and are per- 
fectly fitted by the suggested contact evolution law Eq. ([2]) (k = 1.145\/27r, /i = 
0.55,/? = 0.21). Moreover, if the surface is sufficiently smooth A^ = L/32 we 
do not observe any dependence of results on the Hurst exponent. The evolution 
of the real contact area predicted by Greenwood [18] and BGT [5] theories for 
spectrum bandwidth 7710^4/7/12 = 2 (data from [21]) lie in the confidence interval 
of the results obtained for non-representative surfaces (A; = L) with high Hurst 
exponents H ^ 0.5. Asymptotic limit of Persson theory |23l [23] lies below all 
obtained results. 
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linked to the geometrical represent at ivity, which we define as the proximity 
of the surface heights distribution to a normal distribution [13]. In particular, 
this proximity is important in the range of maximal heights corresponding 
to the only zones which come in contact at small loads. To analyze how the 
representativity depends on the cutoffs, for each pair of cutoff wavenumbers 
kg = L/Xg and ki = L/Xi we generate 30 statistically equivalent surfaces 
using FFT filtering algorithm |2B]. The size of surfaces is L = 1, the dis- 
cretization spacing is AL = 1/1024 and H = 0.8. Fig. [1] represents the 
average L2 error between the heights distribution function of the generated 
surfaces and the corresponding normal distributions G{h, h, a) for different 
cutoffs (see appendix for details). The mean value h and the standard devia- 
tion a = a/ {{h — hy) for the normal distribution G are computed according 
to the generated surface. It follows from this figure that a generated surface 
becomes more and more representative, when the first cutoff wavenumber ki 
increases up to 8 — 16, whereas the variation of kg has a weaker effect on 
the surface representativity if the spectrum is sufficiently rich kg ^ ki. The 
choice Xi = L leads to a surface with a mechanical response which varies con- 
siderably from one surface realization to another. The fact that in average 
this response does not correspond to the response of a bigger rough surface 
with similar random properties is of crucial importance for the study of rough 
contact. The key reasons for this discrepancy are the long-range interactions 
(1/r) between contacting asperities and inevitable boundary conditions on 
lateral sides of the specimen (periodic, symmetric, free). Accordingly, the re- 
sponse of any heterogeneous system with random properties and long-range 
interactions should be studied on a representative system element. 

3 Smoothness of rough surfaces 

Besides the long cutoff wavelength, in generated or experimental surfaces 
there is a limitation connected with their inevitable discreteness, either due 
to a numerical resolution scheme and/or experimental measurements. In real 
self-affine surfaces the power spectrum density decays as a power-law of the 
wavenumber k, C{k) ~ k~'^^^^^\ This law may be preserved down to wave- 
lengths comparable to atomic spacings [2], where the continuum contact me- 
chanics is not valid anymore [27] . To remain in the continuum framework and 
to capture accurately the mechanics, we introduce a short cutoff wavelength 
Xg > AL, which ensures the smoothness of surfaces at a certain magnifica- 
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FIG. 3: Evolution of k with the real contact area fraction A/Aq is depicted for 
different Hurst roughness exponents and cutoff wavelengths A; = L, L/8, L/16, 
a) Xs = L/32, b) = L/128; all results for representative surfaces A; = L/16 are 
fitted by the contact evolution law Eq. ([2]) with k = 1.145\/27r, fi = 0.55, /3 = 
0.21 for (a) A; = L/16, A^ = L/32; with k = 1.145\/2^, /x = 0.56 and /3 € 
[1.120, 1.172] for (b) A; = L/16, A^ = L/128. For comparison, the asymptotic 
limit Kbgt = \/27r from Bush, Gibson, Thomas theory [5] is plotted, while the 
asymptotic limit of Persson theory [231 [23] '^p = \/8/7r 1.60 is out of the plot 
range. 
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tion as, for example, in [7]. The inherent cutoff = AL used, for example, 
in [HI [12] results in a "one node - one asperity" approach, which does not 
allow to correctly reproduce the local change of the contact force with sep- 
aration [13] nor to estimate a realistic growth of localized contact zones but 
only the growth of their number. Considering smooth-rough-surfaces with a 
"truncated self-affinity" ^ AL eliminates these shortcomings. 

4 Description of simulations 

Using the algorithm [26] we generate 12 statistically equivalent surfaces for 
each pair of cutoff wavelengths L/A/ = 1, 2, 8, 16, L/A, = 32, 64, 128, 256, 512 
and different Hurst roughness exponents H = 0.2, 0.36, 0.52, 0.68, 0.84. For 
all the surfaces the constant root mean squared slope VTT^^^I^ = 0.1 is pre- 
served; as before L = 1 and AL = 1/1024. It is important to note that, in 
contrast to geometrical estimations of \J (| V/ip) [6] , which depend on the sur- 
face discretization and thus may underestimate significantly the real value, we 
evaluate the rms slope according to the surface power spectrum [20] . To solve 
the contact problem between a periodic rigid rough surface and a deformable 
fiat half-space we use the spectral based boundary element method [28] with 
some minor improvements [29]. This method allows to solve accurately the 
equations of continuum mechanics under contact constraints and, in contrast 
to asperity-based models, takes into account all underlying mechanics: com- 
plex shapes of asperities, junction of contact zones associated with different 
asperities and long-range deformation of the elastic half-space in response to 
contact forces. For each simulation the contact area fraction reaches ~ 10% 
under the external pressure linearly increasing up to Po/E* = 0.05 within 
30 increments. In Fig. [2] we compare the area-force curves of our numerical 
results and several analytical theories, which were computed and summa- 
rized in [21]. To demonstrate better the nonlinearity of the area evolution 
for surfaces with different cutoffs and Hurst exponents, we present in Fig. [3] 
the proportionality coefficient k expressed from Eq. [T]as k = L'*a/ (|V/iP)/p, 
where p = F/A is the mean contact pressure. 
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5 Results 



The value of k for the observed range (up to 10% of the contact fraction) is 
a decreasing function of the area. For non- representative A; < L/16 or too 
rough \s < 32 AL surfaces, there is a clear tendency (see Fig. [3l b): a higher 
H results in a smaller n. This dependence of the results on the Hurst expo- 
nent corresponds to all up to date numerical investigations of the k constant 
m [161 cm • However, we argue that these results are strongly affected by 
non-representativity and an excessive roughness of exploited surfaces. In con- 
trast to these results, if the mechanics of contact is well resolved (A^ ^ AL) 
and the surface is representative (A; ^ L/16) the real contact area evolution 
does not depend anymore on the Hurst exponent (see dark symbols in Fig. El 
a). This result is in a good agreement with analytical theories and has never 
been obtained before in complete numerical models. Reducing A/ reduces the 
data scatter, which is one of criteria of mechanical representativity. The data 
scatter increases with increasing H, thus, the smoother the surface, the big- 
ger RSSE is needed to capture an average mechanical behavior. The rise in 
K with decreasing A; can be interpreted in terms of Persson theory [231 121] : 
reducing A/ is equivalent to decreasing magnification which rises k. Note 
that our results for k, are not confined between the asymptotic limits of the 
BGT [5] 

'^BGT — V^TT and Persson theories Kp — ^JS/tv [23], as was generally 
observed by other authors [HI [71 [HI [TO] . 

6 Phenomenological Contact Evolution Law 

In light of these results, we propose a new phenomenological contact evolution 
law which is based on our observations of the change in the mean contact 
pressure p = F/A with respect to the applied pressure = F/Aq depicted 
in Fig. [H For considered interval of contact areas A/Aq G [0.01,0.11] the 
decay of dp/dp^ can be approximated by power-law dp/dp^ = f3{AQ/A)'^, 
< /i < 1. Therefore, we obtain (for details see appendix the evolution of 
the real contact area as 



where k is still a proportionality constant between the area and normalized 
force for infinitesimal pressures, two other constants fi and /3 can be easily 
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and uniquely found from the data analysis. For converged results A; = L/\Q, 
As = L/32 by fitting the normalized inverse mean pressure (Fig. [3]) we found 
K = (1.145 ± 0.002)A/27r, which is about 15% higher than the asymptotic 
limit of the BGT theory [5]. From Fig. IHwe uniquely obtain /i = 0.55 ±0.02, 
13 = 0.210 ±0.005. 

7 Conclusions 

To study the growth of the real contact area between rough self-affine man- 
ifolds we obtained statistically meaningful results for different cutoff wave- 
lengths and Hurst exponents. We demonstrated that the real contact area 
evolves nonlinearly with applied force even for reasonably small contact frac- 
tions A/Aq < 0.1. Note that an almost linear evolution was observed for 
rough contact between elasto-plastic materials [211 [I2l [H] , however, in this 
case the Johnson's assumption [17] on replacing two deformable solids by 
one solid with effective elasto-plastic properties and superposed roughnesses 
is not verified and a full simulation of two deformable solids is required. 

In conclusion, we state that to obtain realistic results in rough contact 
analysis, one needs to construct representative surfaces. This representativity 
necessarily requires the long wavelength cutoff to be significantly smaller 
than the specimen size A; ^ L/16. The Hurst exponent does not change 
the evolution of the real contact area with load if both the representativity 
and the smoothness of surfaces are maintained. To describe the evolution 
of the real contact area we proposed a new phenomenological law which 
describes well our results for moderate pressures and reduces to the classical 
BGT law [3] in the limit of infinitesimal contact fractions. However, further 
work is required to verify the universality of this law. In perspective we aim 
to study the pressure distribution and the evolution of the contact area up 
to the full contact, for which the numerical results should be in a better 
agreement with Persson theory [231 121] • An important question, which we 
could not answer in our study, is how the bandwidth parameter influences the 
mechanical behavior of rough surfaces. To address it, one needs to carry out 
numerical simulations on a considerably finer discretization of surfaces than 
was reported here. We hope that our results will motivate new simulations 
based on the notion of representativity and taking into account friction, visco- 
plasticity, adhesion and surface energy. 

The financial support from the European Research Council (ERCstg 
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FIG. 4: The mean pressure derivative with respect to the external pressure 
dp/dpo decreases nonhnearly with increasing fraction of the real contact area 
A/Aq. We plot its evolution for different Hurst roughness exponents and cut- 
off wavelengths A; = L, L/8, L/16 and (a) = L/32, (b) A^ = L/128. For 
representative surfaces the results can be roughly approximated by a power law 
dp/dpo = I3{Aq/A)^; solid line corresponds to parameters /3 = 0.21,^ = 0.55, 
which fit well our results for representative A/ = L/16 smooth surfaces. 
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8 Appendix: representativity of rough sur- 
faces 

The generation of realistic surfaces is the basis for a numerical analysis of 
the contact interaction between solids. It is widely acknowledged that real 
surfaces obey the property of self-affinity [171 [2] over many spatial scales. 
Self-affine surfaces can be described by a power spectrum density (PSD) 
C{k) following a power law of a wavenumber k = L/X 



where H is the so-called Hurst roughness exponent H G (0, 1). In reality the 
PSD may be considered confined between two cutoff wavenumbers: kg and 
ki corresponding to short and long wavelengths and A;, respectively. To 
be representative, a surface has to obey a statistically relevant distribution 
of heights so that, its global mechanical response computed with periodic 
boundary conditions is equivalent to the response of any statistically similar 
rough surface. Since the mechanical response of the generated surfaces is 
affected by their height distributions and the PSD, the surface generation 
procedure has to include both of them into consideration. While the typ- 
ical PSD of self-affine surfaces is easily obtained using a random midpoint 
algorithm [13] or by filtering a white noise [26], the normality of the height 
distribution of simulated surfaces was, to the best of our knowledge, often 
ignored. 

We will demonstrate that even if available numerical techniques permit 
to generate a self-affine surface with a small error on the resulting PSD, the 
distribution of heights strongly depends on the range of wavelengths [A^, A;] 
included in the spectrum. In order to quantify the deviation from a Gaussian 
distribution of heights we introduce the integral L2 error: 



where the mean height h{S), the standard deviation cr{S), and the probability 
density function (PDF) of heights are extracted from surface S and where 



(3) 
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is the number of bins used to evaluate the PDF of the surface. The function 
G stands for the normal distribution 



G{h, h, a) 



exp 



jh-hf 



This measure of the surface representativity has been used in a parametric 
study to identify the effect of the selected wavelength band [A^, A^]. 




FIG. 5: Mean deviation from the Gaussian distribution (computed over 30 
surfaces) with respect to the high cutoff wavenumber ks- Various low cutoff 
wavenumber ki are presented by several curves. The Hurst roughness exponent 
for all surfaces is H = 0.8 and the root mean squared height is rms = 0.1, each 
surface consists of 1024^ points. 



We have synthesized a large set of self-affine rough surfaces (30 per each 
combination of A/ and Ag) containing 1024 x 1024 points by using the filtering 
technique described in [26]. Fig.|5]represents the average error over generated 
surfaces for Hurst exponent H = 0.8 and the number of bins = 500. It can 
be observed that the error decreases with increasing kg for all ki. However, the 
value of ki dictates the order of magnitude of the error for all kg. Therefore 
we may assume that for ki = 16 and kg > 32 the heights distributions of 
generated surfaces do not deviate significantly (in average) from the normal 
one. 
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9 Appendix: contact evolution law 



According to the simulations, the change of the mean contact pressure p = 
F/A with external pressure po = F/Aq decreases approximately as a power 
law of the real contact area fraction A/Aq: 

|^ = /3(^^y, /3>0, 0</i<l (4) 



dpo \ A 

In the limit of small contact we suppose the classical linear relation between 
A/Aq and po 

A K 

However, in contrast to the Bush-Gibson-Thomas model [5] there is no need 
to analyze unreachably small fractions of the real contact area. In the limit 
of the full contact the smoothness of the function is lost 



J ^ = /3, when A ^ Aq and Aq - A > 0, 
\^ = 1, whenA = Ao 



In other words the slope of the real contact area is inversely proportional to 
the contact pressure " needed to establish full contact with the coefficient of 
proportionality d{A/Ao)/dpo = (1 — /3)/Po" however, this theory is supposed 
to be valid only in the range of moderate forces, far from the full contact 
range. Solving Eq. (jlj) gives the following law for the evolution of the real 
contact area 

A _ 1 

A' 



/3 + 



I \ 1— M 



for the considered case we found k = (1.145 ± 0.002)v/27r ^ 2.870 ± 0.005, 
/i = 0.55 ± 0.02, 13 = 0.210 ± 0.005. 

9.1 Details of the CEL derivation 

Starting with the observation that 
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Developing the left part gives 

p = F/A,p^ = F/A^ 



A, 



A-FdA/dF 



dpo A^ 
Substituting this expression in the left part of dHD gives 



A-F^ 

^ ^ dF 

A^ 



/3 Mo"'' 



Ao\A 



(7) 



Grouping terms gives 



F^ = A- PA, 
OF ' ° 



A 



/i-2 



Separating variables and integration of both sides reads as 

A F 



d{A/A,) 



A/Ao-n^/^o) 



f {A/A,Y-' - p 
/i - 1 I {AJA^f-' - /3 



In 



dF 



HF/F,) 



(8) 



Let Fc — 0, — j- and if at very small pressures the real contact area is 
proportional to the normal force with coefficient k, then 



Ar 



^/\VWE* 

Substituting this limit in the bottom boundary of the left integral in Eq. ([8]) 
and integrating, one obtains 

/ \ 



1 



/i — 1 



In 



{A/A, 



(3 



/3 



ln(F/Fe) 



Since /i — 1 < and Fc tends to zero we can neglect /3 in the denominator 



In 



{A/A, 



(3 



■FJA 
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Taking the exponential of both parts one gets: 



{A/A,) 



( 



K 



-Fc/Ao 



Now we can express the real contact area fraction A/Aq as a function of the 
applied pressure po = F/Aq: 
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